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\^ , Aims. We study in some details the statistical properties of the turbulent 2-phase interstellar atomic gas. 

■ Methods. We present high resolution bidimensional numerical simulations of the interstellar atomic hydrogen which describe it 
"^h over 3 to 4 orders of magnitude in spatial scales. 

f») . Results. The simulations produce naturally small scale structures having either large or small column density. It is tempting to 
I ' propose that the former are connected to the tiny small scale structures observed in the ISM. We compute the mass spectrum 
O | of CNM structures and find that Af(M)dM oc M~ 1 ' 7 dM, which is remarkably similar to the mass spectrum inferred for the CO 
^3 clumps. We propose a theoretical explanation based on a formalism inspired from the Press & Schecter (1974) approach and 
^ used the fact that the turbulence within WNM is subsonic. This theory predicts JV(M) oc M" 5/3 in 2D and Af(M) oc M" 16/9 in 
. . 3D. We compute the velocity and the density power-spectra and conclude that, although the latter is rather flat, as observed in 

■ supersonic isothermal simulations, the former follows the Kolmogorov prediction and is dominated by its solenoidal component. 
' This is due to the bistable nature of the flow which produces large density fluctuations even when the rms Mach number (of 
, WNM) is not large. We also find that, whereas the energy at large scales is mainly in the WNM, at smaller scales, it is dominated 

■ by the kinetic energy of the CNM fragments. 

Conclusions. We find that turbulence in a thermally bistable flow like the atomic interstellar hydrogen, is somehow different 
from turbulence in a supersonic isothermal gas. In a companion paper, we compare the numerical results with atomic hydrogen 
observations and show that the simulations well reproduce various observational features. 
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1. Introduction 

The neutral atomic gas (HI) is ubiquitous in the galax- 
ies and is of great importance for the star formation 
process, since molecular clouds form by contraction of 
HI. As such, HI has been extensively observed over the 
years (e.g. Kulkarni & Heiles 1987, Dickey & Lockmann 
1990, Heiles & Troland 2003, 2005, Miville-Deschenes et 
al. 2003). Although HI has also received a lot of attention 
from the theoretical point of view (Field 1965, Field et 
al. 1969, Zel'dovich & Pikel'ner 1969, Pcnston & Brown 
1970, Wolfire et al. 1995, 2003) and in spite of the early 
recognition that HI is a turbulent medium (e.g. Crovisier 
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1981, Heiles & Troland 2005), it is only recently that the 
dynamical properties of HI have been investigated. 

Two lines have been simultaneously pursued by vari- 
ous teams. Attempting to understand the properties of HI 
on large scales, Gazol et al. (2001, 2005), Dib & Burkert 
(2005), de Avillez & Breitschwerdt (2005), have performed 
2D or 3D numerical simulations of turbulent multi-phase 
flows subject to various forcing. On smaller scales, vari- 
ous studies attempting to resolve the physical scales in- 
volved in the problem, have focused on the description of 
the thermal contraction, trying to understand the cooling 
and the condensation of a piece of WNM into CNM ac- 
curately. Hennebelle & Perault (1999, 2000) have studied 
the influence of a converging flow, Koyama & Inutsuka 
(2000) have investigated the influence of a shock propa- 
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gating in a bistable flow whereas Sanchez-Salcedo et al. 
(2002) have considered initially unstable gas. More re- 
cently, Hennebelle & Passot (2006) have investigated the 
influence of Alfven waves propagating into the medium. 
Because of the high numerical resolution required to treat 
the problem, these works were performed in only one 
dimension. Further important extensions of these works 
have been performed in 2 or 3D. Koyama & Inutsuka 
(2002) have generalized their study in 2D and Kritsuk & 
Norman (2002) have computed the evolution of thermally 
unstable gas using 3D simulations. 

In a recent paper (Audit & Hennebelle 2005 here after 
paper I), we have investigated the dynamical evolution of 
colliding flows of WNM, paying particular attention to the 
effect of turbulence. The main conclusions of this study 
are as follows: first, we confirm, in a context where the 
flow is strongly turbulent, the result obtained by Koyama 
& Inutsuka (2002) in their study of shock propagation 
in WNM. The CNM is very fragmented in small clouds 
which have a dispersion velocity equal to a fraction of 
the WNM sound speed. Second, as in Gazol et al. (2001), 
although at much smaller scales, we find large fractions 
of thermally unstable gas and we show that turbulence 
is able to stabilize transiently the thermally unstable gas. 
Third, we show that in spite of the strong turbulence, the 
CNM fragments are generally not very far from pressure 
equilibrium with the surrounding WNM. Colliding flows 
of atomic interstellar gas have been further investigated 
by Heitsch et al. (2005, 2006) and by Vazquez-Semadeni 
et al. (2006). 

One of the crucial issues in modeling atomic hydrogen 
is the numerical resolution which is very demanding (see 
section |2~3|) . In this paper, we present the result of very 
high resolution numerical simulations allowing description 
of scales ten times smaller than in paper I. We perform 
a set of simulations with various numerical resolutions to 
study its influence. We also study the influence of thermal 
conduction by varying its value. The high resolution al- 
lows a better description of the flow down to much smaller 
scales and offers the possibility to perform various statis- 
tical studies of the flow and CNM structures properties. 

In a companion paper (Hennebelle ct al. 2006, paper 
III), we make preliminary comparisons between the 2D 
simulations presented here and various observational re- 
sults. In particular, the high resolution permits the com- 
parison with various small scale structures that have been 
observed in HI, like the so-called tiny small atomic struc- 
ture (TSAS) and recent low column density clouds ob- 
served by Braun & Kanekar (2005) and Stanimirovic & 
Heiles (2005). 

In Sect. 2, we describe the numerical setup and method 
and discuss briefly the various scales which have to be 
adequately treated. In Sect. 3, we present one typical 
timcstcp of the simulation and discuss qualitatively var- 
ious features of the flow which are worth quantifying. 
Section 4 presents statistical density and the pressure 
PDF. Section 5 presents density and velocity power- 
spectra as well as the energy spectrum. In Sect. 6, the 



mass spectrum of the CNM structures in the simulation 
is computed and in Sect. 7, a theoretical explanation is 
proposed. Section 8 concludes the paper. 

2. Initial conditions and method 

The initial conditions and the methods are very similar 
to those of paper I. Here, we describe them briefly for 
self-consistency. 

2.1. Equations, notations and numerical method 

We consider the usual fluid equations for a radiatively 
cooling gas including thermal conductivity namely, 

d t p +v.H =o, (i) 

d t pu + V.[pu®u + P]=0, (2) 
d t E + W.[u{E + P)} =-£(p,T)+V(k(T)VT). (3) 

p is the mass density, u the velocity, P the pressure, E 
the total energy and C the cooling function which in- 
cludes Lyman-Q, C + and O line cooling and grains pho- 
toelectric heating. The gas is assumed to be a perfect 
gas with 7 = 5/3 and with a mean molecular weight 
p = IAttih, where nin is the mass of the proton, k is 
the thermal conductivity and is given by k{T) = "fC v r](T) 
where C v = k h /m H /{l - 1), V = 5-7 x 10" 5 (T/l K) 1 / 2 g 
cm -1 s _1 and kb is the Boltzmann constant. 

We use the HERACLES code to perform the simula- 
tion. This is a second order Godunov type hydrodynamical 
code. Godunov methods are now widely used and tested 
and are particularly well suited to the treatment of shocks 
(e. g. Toro 1997). The size of the computational domain is 
20 pc and the resolution ranges from 600 2 to 10000 2 cells 
leading to a spatial resolution ranging from 3.3 x 10~ 2 pc 
to 2 x 10~ 3 pc. 

2.2. Initial and boundary conditions 

The boundary conditions consist in an imposed converging 
flow at the left and right faces of amplitude 1.5 x C StWnm 
on top of which turbulent fluctuations of amplitude e = 2 
(see paper I for accurate definitions) have been superim- 
posed. On top and bottom faces, outflow conditions have 
been setup. This means that the flow is free to escape the 
computational box across these 2 faces. 

The initial conditions consist of uniform WNM (n ~ 
0.8 cm -3 ) at thermal equilibrium. It is worth stressing 
that initially no CNM is present in the computational box. 
The simulations are then runned until a statistically sta- 
tionary state is reached. Typically, this requires about 5 
to 10 box crossing times. At this stage, the incoming flow 
compensates on average the outgoing material. 

For the highest resolution runs, this would lead to very 
long simulations and we therefore adopt a slightly different 
strategy. We start the run with a lower resolution of 1250 2 
cells and wait until the statistically stationary regime is 
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Fig. 1. Density and velocity fields at one time step of the 10000 2 cells simulation corresponding to a spatial resolution 
of 2 x 1CP 3 pc. Note that because of the large dynamics, the arrows are not easily seen in the CNM clouds. 



obtained. Then, we increase the resolution of the simu- 
lation by a factor 2 and restart the now 4 times larger 
simulation. This new simulation is runned until the sta- 
tionary regime is reached. Typically this takes about two 
box crossing time. This procedure is then repeated until 
the desired spatial resolution is obtained. 

2.3. Scales 

Contrary to polytropic fluids, a 2-phase medium has vari- 
ous spatial and temporal characteristic scales which have 
to be adequately described. Here, we briefly recall them. 



2.3.1. A static spatial scale 

The Field length (Field 1965) is the length at which ther- 
mal diffusivity becomes comparable to the heating and 
cooling term. Its typical value is about 0.1 pc in the WNM 



and about 10~ 3 pc in the CNM. The Field length is also 
the typical scale of the thermal fronts that connects the 
cold and warm phases. As demonstrated by Field (1965), 
it determines the smallest wavelength at which thermal 
instability can grow and, therefore, the size of the small- 
est CNM structures. The effect of varying the thermal 
diffusivity is discussed in the following. 



2.3.2. Three dynamical scales 

We distinguish three scales whose origin is due to dynam- 
ical processes. The first scale is the cooling length of the 
WNM, i.e. the product of the cooling time and the sound 
speed within WNM, A coo i = t coo i x C s .wnm- This scale 
is about 10 pc and corresponds to the typical length at 
which WNM is non linearly unstable and can be dynami- 
cally triggered by a compression into the unstable regime. 
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Since the ratio between the CNM and WNM den- 
sity is about hundred, the size of a CNM structure 
formed through the contraction of a piece of WNM of size 
Acoob will be typically, hundred time smaller (assuming 
a monodimcnsional compression). Therefore, the size of 
the CNM structures is about the cooling length of WNM 
divided by hundred, leading to about A coo i/100 ~0.1 pc0. 

Finally, as first pointed out by Koyama & Inutsuka 
(2002), the fragments of CNM have a velocity dispersion 
equal to a fraction of the sound speed of the medium in 
which they are embedded, i.e. WNM. Since the contrast 
between the sound speed within the two phases is about 
10, it means that CNM structures undergo collisions at 
Mach number, M ~ 10. If, for simplicity, we assume 
isothermality and apply Rankine-Hugoniot relations, we 



1 In fact, as we will see in the following, the size of the CNM 
structures follows a density distribution. The numbers given 
here are simply indicative 



obtain that the size of the shocked CNM structures is 
given by the size of the CNM structures divided by M 2 
which is about 10 -3 pc, whereas the density peak is given 
by 100 x M 2 ~ 10 4 cm~ 3 . Note that since the effective 
polytropic index of the CNM, J e ff, is closer to 0.7 than to 
1, these numbers are underestimated by a factor of about 
7. 

The consequences of not solving properly these scales 
have been previously discussed qualitatively in paper I. In 
the following, quantitative estimates are given. 

2.4. Comparison with other simulations 

As explained in the previous section, it is actually not pos- 
sible to treat simultaneously all the scales involved in the 
physics of the atomic gas. This implies compromises and 
induces various choices of initial, boundary and forcing 
conditions. In order to make clear the domain of valid- 
ity of various related works, we now discuss the numerical 
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Fig. 3. Logarithm of the pressure field (in K cm' 3 ) corresponding to Fig. [2] 



and physical setups which have been chosen. It is worth 
stressing that the different choices are complementary and 
equally interesting since they allow one to focus on differ- 
ent aspects of the same physical problem. 

Koyama & Inutsuka (2002) have investigated the prop- 
agation of a shock into WNM. They have an extremely 
high spatial resolution (7 x 10~ 4 pc) and a relatively mod- 
est box size of about 0.25 pc (in the transverse direction). 
With these choices, the CNM structures are very well re- 
solved but the statistics are low and the structures are 
small. In a recent paper (Koyama & Inutsuka 2006), they 
consider a larger box and try to understand the proper- 
ties of the turbulence that develops without any external 
forcing. 

Hcitsch et al. (2005, 2006) and Vazquez-Semadeni et 
al. (2006) have focused on the formation of a CNM turbu- 
lent layer. They consider a simulation box of about ~10 pc 
and have a resolution of about 10 -2 pc. Such a setup which 
is very similar to the one used in paper I, allows to treat 



reasonably well the large scales, but does not provide a de- 
scription of the CNM structures as accurate as the CNM 
structures obtained in Koyama & Inutsuka (2002). Also 
in these work, the authors focus on the transient regime, 
rather than on a time-independent situation. 



Gazol et al. (2001, 2005) and Piontek & Ostriker (2004) 
have focused on larger simulation boxes (100 pc-1 kpc) 
and coarser spatial resolution (0.5-5 pc). The same kind of 
scales has been investigated by De avillez & Breitschwerdt 
(2005) using an AMR scheme. This leads to a much bet- 
ter description of the larger scales and a more realistic 
self-consistent forcing (e.g. due to supernovae or the de- 
velopment of the magneto-rotational instability) that the 
choice makes in this work, but only marginal description 
of CNM structures can be achieved. 
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3. Results 

In this section, we present results for the 10000 2 cells sim- 
ulation which corresponds to the highest numerical reso- 
lution runs. We discuss qualitatively the main properties 
of the 2-phase turbulent flow and introduce the most im- 
portant aspects which are studied quantitatively in the 
following and in paper III. 

3.1. Main features 

Figure Q] displays the density and velocity fields of the 
simulation at one representative timestep (i.e. well after 
statistically stationary state has been reached) whereas 
Fig. [U shows a spatial zoom of 5 pc. 

With our choice of boundary conditions which entails 
a converging flow, the formation of a layer-like structure 
is triggered in the middle of the computational box (x ~ 5 
pc, y ^ 2 - 18 pc). 

The fragmentation of such a layer due to the develop- 
ment of various instabilities, like the Vichniac instability 
(Vichniac 1994) or the Kelvin-Helmholtz instability has 
been studied by Walder & Folini (1996, 1998) and Heitsch 
et al. (2005, 2006) for a radiatively cooling gas and by 
Folini & Walder (2006) for strongly supersonic isother- 
mal gas. Here, since the converging flow imposed at the 
boundary presents strong velocity fluctuations, the tur- 
bulence and the fragmentation of the layer are not easily 
attributable to a single instability and at least, partly due 
to the driven turbulence. Indeed, in paper I, we vary the 
amplitude of the shear velocity fluctuations which are su- 
perimposed to the converging flow. As expected, we find 
that the stronger these fluctuations, the stronger the tur- 
bulence in the computational domain. Near the bound- 
aries of the computational box, the structure of the flow 
is somehow different because the converging flow is con- 
centrated in the middle of the computational box faces, 
whereas the ram pressure at the edge of the box faces 
almost vanishes. 

The overall structure of the flow appears to be highly 
complex. The cold phase is very fragmented and the two 
phases are highly intervowen. There are large regions of 
space where there is only WNM and regions where most 
of the mass is obviously in the CNM. However, even in 
these denser regions, WNM is still present and intrusive. 
Only few large CNM structures have formed, most of them 
appear to be very small. This is more clearly seen in Fig. [5] 
which also reveals that the structures are often connected 
to each other by filaments of lower density and that the 
density of the CNM varies over 2 orders of magnitudes. 

3.2. Shocked CNM: TSAS ? 

While most of the CNM structures have a density of about 
100 cm -3 (e.g. structure located at (x,y) = (3,8) pc), 
shocked regions present densities up to ~ 10 4 cm~ 3 (e.g. 
structure located at (x,y) — (1.66, 7.07) pc). This is more 
clearly seen in Fig. |4]which shows a spatial zoom of Fig. [2] 



and in Fig. [5] where a cut along the y-axis is displayed. 
These two figures also reveal that the density peak corre- 
sponds to a peak of pressure reaching a value of about 10 5 
K cm~ 3 and that this strong fluctuation has been induced 
by a converging flow of velocity ±5 km/s. Altogether, this 
is perfectly consistent with the orders of magnitude given 
in paper I and recalled in Sect. 12.31 to characterize the 
shocked CNM structures. The density (~ 10 4 cm -3 ) and 
the size (~ 400 AU) appear to be close to what is inferred 
for the TSAS structures observed in the CNM (e.g. Heiles 
1997), although the more extreme events present densi- 
ties up to 10 times this value and scale roughly 10 times 
smaller. It is therefore tempting to propose that TSAS 
are produced by supersonic collisions between CNM frag- 
ments. It should be noted that the density peak is not 
fully resolved because of insufficient numerical resolution 
and therefore it is expected that this event should have 
indeed a larger peak density and a smaller size if the nu- 
merical resolution was higher. It is also important to em- 
phasize that the cooling function used in this simulation 
is probably not accurate for such high densities. 

3.3. Small column density clouds 

In a similar way, the properties of the smallest structures 
that form in the simulations appear to be reminiscent 
of the tiny structures which have been recently observed 
by Braun & Kanekar (2005) and Stanimirovic & Heiles 
(2005). Many of them can be seen in Fig. [2] (e.g. for y < 6 
pc) and a cut through two small CNM structures is dis- 
played in Fig. [5] (y ~ 6.68 and 6.83 pc). The density of 
these two structures is about 30 and 100 cm" 3 whereas 
their size is about 10~ 2 pc. 

Further examples are shown in paper III in which syn- 
thetic HI spectra are also discussed. 

3.4. Pressure distribution 

Another important aspect is revealed by the pressure field 
displayed in Fig. [31 The pressure field appears well orga- 
nized with large scale systematic gradients. At large scales, 
the pressure and the density fields appear to be broadly 
anti-correlated. Indeed, the central cloud (located roughly 
between x = 1 and x — 4 and y — 6 and y = 9) is 
surrounded by high pressure WNM, whereas its internal 
pressure is a few times lower. This is a consequence of the 
radiative cooling. When the WNM fluid particles enter 
into the shocked layer, their internal pressure goes up like 
in an adiabatic shock, then they cool and finally condense 
out. Indeed, the cooling length of shocked (Mach 1.5 — 2) 
WNM is few parsec which is the length of the high pressure 
WNM layer seen on the left part of the box. Interestingly, 
the high density shocked CNM is located at the interface 
between the high pressure WNM and the low pressure 
cloud. Note that since the WNM pressure (P ~ 10 4 K 
cm~ 3 ) is lower by at least one order of magnitude than 
the large pressure fluctuation (P ~ 10 5 K cm -3 ), it can- 
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not be the triggering agent. However, this high pressure 
WNM is in part responsible of these large pressure fluctu- 
ations since it actively induces the formation of relatively 
high velocity CNM structures, which then undergo super- 
sonic collisions. 

At small scales, the situation is much different. The 
pressure is poorly correlated with the density field dis- 
played in Fig. [2 Most of the CNM structures do not have 
high thermal pressure (see Fig. [5]) with the notable excep- 
tion of the densest regions which correspond to shocked 
CNM. Indeed, some of the CNM structures have a thermal 
pressure lower than the surrounding WNM as for example, 
the 2 structures located at y = 7.15 and 7.3 pc in Fig. \5\ 
This low pressure indicates that the structures have just 
formed, so that there was not enough time for the initially 
low thermal pressure induced by the radiative cooling to 
readjust to the surrounding higher pressure. 

Finally, the density field indicates that the CNM struc- 
tures are well identified and connected to the surrounding 



WNM by very stiff fronts. Moreover, as in the static 2- 
phase medium (Field et al. 1969, Wolfire et al. 1995), most 
of the CNM structures are not far from quasi-pressure 
equilibrium with the surrounding medium (in the sense 
that pressure fluctuations are small with respect to the 
density contrast between the two phases). This is likely 
a consequence of the large density contrast between the 
phases. First, with a density contrast of about hundred, 
the sonic waves propagating in the WNM as well as the 
eddies get mostly reflected when they reach a WNM/CNM 
interface, probably making the energy injection inside the 
structures relatively inefficient. Second, since the scale of 
the CNM structures is small compared to the scale of the 
WNM flow and since most of the turbulent energy is on 
the largest scales, the WNM flow is relatively uniform at 
the scale of the CNM structures and mostly advects them 
uniformly. 

To summarize, it seems that turbulent 2-phase flows 
cannot be accurately described as a polytropic turbulent 
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Fig. 5. Density, pressure and velocity profiles at x = 1.65 
pc for y ranging from 6.6 to 7.6 pc. Note that this corre- 
sponds to the line of sight which intercepts the cell having 
the highest density at this timestep. 

flow, neither as a static 2-phase medium, but, instead, 
they present a kind of duality, some aspects being more 
reminiscent of turbulent flows and other of static 2-phase 
media with fluctuations which can be up to hundred times 
the mean CNM density or pressure. In the following, we 
therefore quantify the mass distribution of CNM struc- 
tures, as well as the density and pressure distribution of 
the flow. We also study the statistics of the flow by com- 
puting various power-spectra, in order to characterize its 
scale dependence. In paper III, we quantify further the 
numerical results by computing line of sight statistics and 
various properties of the CNM structures. 

4. Density and pressure probability distribution 
function 

Figures [5] and [7] show the probability distribution function 
of the logarithm of the density and pressure respectively. 
Full black lines correspond to the runs having I 2 = 10000 2 
cells, whereas the dotted red lines are for the runs with 
lg = 5000 2 cells. The density pdf presents two peaks which 
correspond to the WNM and CNM densities. The pdf val- 
ues at densities between 5 and 50 indicates that there is 
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Fig. 6. Probability distribution function of the logarithm 
of the density field. 
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Fig. 7. Probability distribution function of the logarithm 
of the pressure field. 

a significant fraction of thermally unstable g first 
stated by Gazol et al. (2001) and confirmed in paper I. 

In principle, the density pdf allows estimation of the 
fraction of the dense gas that has been compressed by 
the ram pressure. However it is seen that the results for 
the high density part do clearly depend on the numerical 
resolution, meaning that higher resolution runs will lead 
to a higher fraction of high density gas. It should also be 
clear that since the fraction of high density gas depends 
on the cooling processes, it will certainly be modified by 
a more accurate treatment of the microphysics. Having 
these limitations in mind, we can indicatively say that the 
fraction of gas denser than ~ 10 3 cm" 3 is about 1-3 % 
for this simulation. It is important to keep in mind that, 
as emphasized in paper I and in Gazol et al. (2005), this 
obviously depends on the driving of the turbulence. 

The pressure pdf indicates that the average pressure is 
about 6000 K cm" 3 , which is the pressure of the WNM in- 
jected at the boundary. The large pressures are induced by 
the ram pressure of the incoming flow. The pressure dis- 
tribution is reminiscent of pdf found by various authors 
in different contexts (Passot & Vazquez- Semadeni 1998, 
Scalo et al. 1998). Comparison with the results presented 
by Gazol et al. (2005) is more straightforward since they 
consider a thermally bistable flow as well. As in their case, 
we find that a near-powerlaw tail develops at high and low 
pressure (in their case they report that this depends on 
the driving length). Because of the higher numerical res- 
olution, we reach higher values of pressure. Our pressure 
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Fig. 9. Compensated power-spectrum of velocity (full 
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line), compressible component (dashed line). The spectra 
have been multiplied by k a where the value a is indicated 
on the figure. 
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Fig. 12. Energy spectrum for numerical resolutions equal 
to dx = 8 x 1(T 3 , 4 x 1(T 3 and 2 x 10 -3 pc. 




Fig. 13. Power-spectrum of the square root of density field 
for data obtained by taking min(p, po) for various values 
of p . 



Fig. 10. Isothermal simulation. Compensated power- 
spectrum of velocity (full line), solenoidal or shear com- 
ponent of velocity (dotted line), compressible component 
(dashed line). The spectra have been multiplied by k a 
where the value a is indicated on the figure. 



pdf appears to be somehow similar to the results presented 
in their figure 11 corresponding to a large scale driving. 
However, the significant difference between the 2 numeri- 
cal experiments precludes detailed comparisons. 
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Fig. 14. Energy spectrum for data obtained by taking 
min(p, po) for various values of Pq. 
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5. Velocity and density power-spectra and energy 
spectrum 

5.1. Incompressible and isothermal turbulence 

An important diagnostic commonly used in fluid dynam- 
ics to characterize the scale-dependent structure of the 
flow is the power-spectrum of the fluid variables such as 
the velocity and the density fields. The power-spectrum 
equal to the square of the module of the Fourier trans- 
form, depends on all components of the wave vector k. 
The shell average power-spectrum is the averaged of the 
power-spectrum in spherical shells in k-space. In the fol- 
lowing, the power-spectrum will refer to shell averaged 
power-spectrum. 

In incompressible turbulence, the velocity power- 
spectrum, Py, is related to the energy spectrum as 
E(k) oc k D ~ 1 Py, where D is the space dimension. The 
Kolmogorov (Kolmogorov 1941) prediction for the energy 
spectrum of homogeneous, 3D and incompressible turbu- 
lence is: E(k) oc fc~ 5 / 3 . Note that this corresponds to 
power-spectrum proportional to Py oc fc~ n / 3 for 3D flows 
and to Py oc fc -8 / 3 for 2D flows. However, the behaviour of 
2D incompressible turbulence is very different from the 3D 
turbulence because of enstrophy conservation (Kraichnan 
1967). At scales larger than the injection scale, one expects 
to recover the Kolmogorov energy spectrum whereas at 
scales smaller than the injection scale, Kraichnan (1967) 
predicts kPy oc fc -3 . 

In ID compressible turbulence, when shocks dominate, 
a Burgers-like spectrum Py oc k~ 2 is expected (Elsasser & 
Schamel 1976). Passot et al. (1988) have confirmed, using 
isothermal bidimensional numerical simulations, that the 
compressible part of the velocity field follows the same be- 
haviour whereas for the solenoidal part, the prediction of 
Kraichnan (1967) is well verified, kPy oc k~ 3 . Kim & Ryu 
(2005) recently explored numerically in some details the 
3D case. They show that the density power-spectrum (P p ) 
for low rms Mach number flow is close to the Kolmogorov 
power-spectrum, whereas it is much flatter at high rms 
Mach number, i.e. k 2 P p oc k~ a with a ~ 0.5 for an rms 
Mach equal to 12. They interpret this result as the con- 
sequence of the density distribution being highly concen- 
trated in sheets and filaments. Flat density power-spectra 
have also been observed in MHD isothermal simulations. 
Padoan et al. (2004) report a = 0.25 for a model with 
equipartition between magnetic and kinetic energy and 
a = 0.75 for a super- Alfvenic model (see also Beresnyak 
et al. 2005). 

5.2. Velocity power-spectrum 

Figure [8] shows the velocity power-spectrum for various 
numerical resolutions, whereas Fig. [9] shows the compen- 
sated power-spectra (i.e. k a Py the value of a being given 
in the figure) of the solenoidal and compressible compo- 
nents of the velocity field. For reference the Kolmogorov 
(1941) and the Kraichnan (1967) prediction have been dis- 



played in Fig.[H As in Passot et al. (1988), the compress- 
ible part has a power-spectrum kPy oc fc~ 2 . The veloc- 
ity power-spectrum, as well as the power-spectrum of the 
solenoidal component, appear to be in good agreement 
with the Kolmogorov law, kPy oc fc~ 5 / 3 for k between 
~-3.3 and ~-1.8. 

This may be at first surprising because, as recalled 
above, the prediction of Kraichnan for 2D incompressible 
fluid is that there is an invert cascade above the energy in- 
jection scale, leading to the velocity power-spectrum pre- 
dicted by Kolmogorov and a direct enstrophy cascade at 
scale below the injection scale leading to kPy oc fc~ 3 . 
However, this result is a consequence of the conservation 
of vorticity in incompressible or barotropic flow. When the 
flow is non barotropic the vorticity equation writes: 

d t 0J + V x (u x v) = VP x V(-) (4) 

P 

In a radiatively cooling flow, as the interstellar atomic hy- 
drogen, the pressure and the density are not necessarily 
correlated as already emphasized in section 3. This means 
that the baroclinic term on the right hand side, is not 
small with respect to the second term of the left hand 
side, implying that enstrophy is not conserved. To verify 
this, we have computed the mean value of the (norm of 
the) baroclinic term and of the second lhs term over the 
whole simulation. We find that the former is on average 10 
times smaller than the later. We also find that for about 
7% of our grid points, the ratio between the two is higher 
than 0.5. It seems likely to us that this constitutes a suffi- 
ciently large deviation from the strict enstrophy conserva- 
tion to explain the deviation from Kraichnan's prediction. 
Unsurprisingly, we find that the higher values of the baro- 
clinic term are obtained within the cold structures where 
the density changes much but not the pressure. 

To confirm the influence of the baroclinic term, we 
have performed isothermal simulations with exactly the 
same setup as simulations with cooling (same code and 
same initial and boundary conditions). In these simula- 
tions, the baroclinic term vanishes exactly. The compen- 
sated spectra are shown in Fig. [TDJ The power-spectrum 
of the solenoidal component is very stiff and roughly pro- 
portional to k~ whereas the compressible component 
has a power-spectrum proportional to about k . These 
numbers are indeed much closer to the results of Passot 
et al. (1988) and to the prediction of Kraichnan (1967). 

Finally, we note that, as revealed by Fig. [9l although 
compressible and shear modes are comparable at large 
scales, the second dominates at small scales for log(fc) > 
—3 (solid and dotted lines are very close between — 1 > 
log(fc) > -3). 

5.3. Density power-spectrum 

Figure [IT] shows the power-spectrum of the density field. 
It appears to be flat since one finds kP p oc k~° . This 
is reminiscent of what is observed by Kim & Ryu (2005) 
although in the present case the strong and stiff density 
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fluctuations are induced by the 2-phase nature of the flow 
rather than by highly supersonic motions. 

Kim & Ryu (2005) proposed that the flat density 
power-spectrum can be understood by the fact that the 
Fourier transform of a Dirac function is independent of 
k leading to P p oc fc° in ID. Wc believe that this expla- 
nation is also valid in the present case and that the den- 
sity power-spectrum is related to the size distribution of 
the structures (Fig. [T5lfT7|) . Indeed since the structures are 
bounded by contact discontinuities, they have stiff bound- 
aries and they resemble, depending on the scales at which 
one is considering them, either to the Dirac or Heaviside 
functions. 

As a matter of fact, the end of the inertial range which 
occurs at about k ~ 0.1 (about 10 cells in the simulation) 
corresponds to the length at which the structure distri- 
bution is influenced by the numerical diffusivity implying 
that the number of density peaks at these scales becomes 
gradually smaller. 

Interestingly, Desphande et al. (2000) using 21cm line 
absorption data probing linear scales between 10 -2 and 
3 pc, therefore comparable to the present simulation, re- 
port density power-spectrum k 2 P p oc k~ a , in two differ- 
ent regions for which a = 0.75 and a = 0.5. These values 
are significantly flatter than the Kolmogorov index and at 
least for the second one, are very close to the value of 0.4 
which is inferred in our simulations. 

5.4. Energy spectrum 

Figure fl2l shows the power-spectrum, Pe, of the quantity 
yfpv. The quantity fc-Pg corresponds to the energy spec- 
trum in the incompressible case and following Klcsscn ct 
al. (2000), we call it the energy spectrum, E{k). It is rather 
different from the Kolmogorov prediction and significantly 
flatter. In the inertial range E{k) oc fc -0 - 9 , meaning that 
the energy is approximately equally distributed in the k- 
space. This appears to be a direct consequence of the den- 
sity power-spectrum discussed above which relies on the 
strong density fluctuations due to the 2-phase nature of 
the flow. In order to confirm this, we have calculated the 
energy spectrum and the power-spectrum of the square 
root of p for data obtained by taking min(p, po)- The re- 
sults for various values of po, namely 2, 3 and 10 cm~ 3 as 
well as the result for the untruncated field, are shown in 
Figs. [13] and [14] When po decreases, P becomes stiffer 
(note that for the untruncated field, we find kP ^ oc k~ - 6 
in the inertial range) . The same behaviour is also observed 
for the energy spectrum which becomes closer to the ve- 
locity power-spectrum shown in Fig. [5) 

For k larger than ~ 5x 10~ 3 (corresponding to a scale 
of about 1 pc) and po < 3 cm" 3 , the energy spectrum 
becomes significantly lower than the energy spectrum ob- 
tained for the untruncated field. This clearly indicates that 
at scales smaller than about 0.5 pc, most of the energy 
is stored in the CNM bulk motions. In other terms, in 
our numerical experiment, the 2-phase nature of the flow 
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Fig. 15. Number of structures per logarithmic interval of 
mass (in solar mass per parsec). Various numerical reso- 
lutions are shown. 
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Fig. 16. Total mass within CNM structures per logarith- 
mic interval of mass. 

appears to play an important role and constitutes a sig- 
nificant deviation of single phase hydrodynamics at scales 
smaller or comparable to about 1 pc. 

6. Mass spectrum 

As in paper I, we extract the structures by a simple clip- 
ping algorithm. Since the cold structures are connected 
to the surrounding warm gas by stiff thermal fronts such 
simple criteria has a clear physical meaning. The threshold 
is equal to 30 cm -3 which lies in the thermally unstable 
regime (note that we have verified that taking 10 cm~ 3 
instead does not change the results significantly). 

6.1. Numerical result 

Figure [TBI shows the number of structures, n s , per log- 
arithmic mass interval, dn s /d\ogM (note that since the 
simulations are 2D, the mass is expressed in solar mass 
per parsec). The largest structures have a mass of about 
one solar mass per parsec whereas the mass of the small- 
est structures is about ~ 10~ 5 solar mass per parsec. The 
number of structures, M{M), decreases with mass, and fol- 
lows approximately: dn s / dlogM oc 1/A-f 3-1 with [3 ~ 1.7, 
until log(Af) ~ —3. This implies that the number of struc- 
tures per mass interval is dn s /dM = Af(M) oc 1/M 13 . 

For log(M) < —3, dn s /d\ogM increases with mass. 
This part of the distribution is most likely affected by 
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Fig. 17. Total number of cells per logarithmic interval of 
cells number (expressed in cells corresponding to the high- 
est resolution run). 

the numerical resolution. This is confirmed by the fact 
that the numbers of low mass objects decreases with the 
numerical resolution. Interestingly enough, the number of 
massive structures becomes independent of the resolution 
for dx < pc (see solid, dotted and dashed curves). 

This indicates that numerical convergence is reached for 
the large structures (log(Af) <; —3). It is however clear 
that no numerical convergence has yet been reached for 
the smallest structures. We also note that the number of 
structures more massive than ~ 0.3 M s pc -1 , appears to 
be smaller than the value expected from the power-law 
fitted above. This is most likely due to the finite size of 
the numerical experiment. 

For log(M) > —3, the total mass contained in 
a logarithmic interval of mass, follows approximately 
dM tot /d\og oc M 2 Af ~ M 2 -P ~ A/o.25-0.3 which indi _ 
cates that, although the mass of CNM is, in principle, 
dominated by the large clouds, the index of the power- 
law is nevertheless very shallow. This is confirmed by 
Fig. [TBI which shows dM to t/d\og as a function of log(M). 
It is seen that the mass in large structures is only slightly 
larger than the mass in hundred to thousand times smaller 
structures. This indicates that the description of the small 
structures is very important in order to achieve a fair de- 
scription of the flow and its dynamics. 

Complementary information is given in Fig. [17] which 
shows the total number of cells, dn c ^ ot / d\ogn Cl per loga- 
rithmic interval of computational cells number, n c , con- 
tained in the structure (expressed in cells of the highest 
resolution run). The physical size of the structures is ob- 
tained by taking the square root of the cell number and 
multiplying by 2 x 1CP 3 pc which is the numerical resolu- 
tion of the 10000 2 cells run. The size at which the structure 
number starts to decrease because of insufficient numerical 
resolution is about 10 (~ ^frT c for n c ~ 100). The length 
of the bigger structures is roughly 300 cells (~ ^Jn~ c for 
n c ~ 10 5 ) which corresponds to about 0.6 pc. 

6.2. Influence of thermal conduction 

The importance of thermal conductivity has been recently 
pointed out by Koyama & Inutsuka (2004). They show 



1 : 




-4 EL 



-6 -4 -2 

log(M) (M,/pc) 

Fig. 18. Same as Fig. [TB] except that the influence of var- 
ious conductivity is shown. 

explicit cases for which the number of structures which 
forms do depend on the conductivity. In order to inves- 
tigate the exact influence of the thermal conductivity on 
the structures, we have performed various runs in which 
its value has been changed. The results are presented in 
Fig. [HI Four cases have been explored. The full black line 
is for standard ISM conduction, the dashed green line is 
when no explicit conduction is taken into account and the 
dotted blue line is for a conduction 10 times the standard 
ISM conditions. Finally, in order to determine the impor- 
tance of the dependence in temperature of the conduction, 
we also perform a run in which the conduction does not 
vary with temperature and is equal to the conduction at 
T = 8000 K. The corresponding run is shown by the dot- 
dashed red line. The number of structures more massive 
than ~ 10~ 3 solar mass per parsec appears to be indepen- 
dent of thermal conduction. There is a clear influence of 
thermal conduction on smaller structures, larger conduc- 
tion runs produce less small structures than smaller con- 
duction runs. This is qualitatively in good agreement with 
the trend expected from Field's linear analysis that shows 
that the conductivity prevents the formation of structures 
smaller than the Field length. However, the discrepancy 
between the different curves is not very large. We therefore 
conclude that thermal conduction, although not insignifi- 
cant, does not appear to play a dominant role. 

7. Theoretical derivation of the mass spectrum 

Cloud mass spectrum is an important feature of the flow. 
It would be useful to have a theoretical model to explain 
the numerical results. Here, we propose an approach in- 
spired from the Press & Schecter (1974) formalism (see 
also Padmanabhan 1993). 

7.1. Principle and justification of the analysis 

The starting point of the Press & schecter (1974) formal- 
ism is to consider the density spectrum of the initial fluc- 
tuation that eventually lead to galaxies. They consider the 
density field smoothed at a given scale, R, by a window 
function W and assume that the probability of finding a 
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Fig. 19. Isothermal simulation. Density field smoothed at 
scale R — b x dx, for various values of b. Solid line is 
b = 100, dashed line is b — 200, dashed line is b — 500, 
whereas dot-dashed line is b — 1000. 
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Fig. 20. 2-phasc simulation. Density field smoothed at 
scale R = b x dx, for various values of b. Solid line is 
b = 100, dashed line is b — 200, dashed line is b — 500, 
whereas dot-dashed line is b — 1000. 

contrast density 5 = p/p — 1, p being the mean density, is 
where 

<r 2 {R) = Js 2 (k)W%{R)d 3 k, (6) 

(see e.g. Padmanabhan 1993, chapter 5). In this last ex- 
pression 8 2 {k) is the power spectrum of S. 

Therefore, the first question that must be addressed 
is whether this assumption is justified in the case of the 
atomic hydrogen. More precisely, since the CNM struc- 
tures, are the final product of density fluctuations aris- 
ing in the WNM, the question which has to be ad- 
dressed is whether the spectrum of the density fluctua- 
tions smoothed at a scale R, is reasonably described by 
the probability stated by Eq. [5] and what is the value of 
<t(R). One difficulty to address this question is that in the 
multiphase simulations presented previously, the WNM 
and the CNM are mixed together. They are therefore not 
representative of the initial state of the WNM which leads 
to the CNM structures. On the other hand, the early times 
of the simulation reflect the numerical setup rather than 
the initial conditions. To overcome this difficulty, we con- 
sider the simulation with cooling presented above but in 



which we have removed the CNM by taking min(p,po), 
with po = 3 cm~ 3 (see section 5.4). As shown by Fig. [T3l 
the power-spectrum of ^fp is nearly Kolmogorov. We have 
checked that the power-spectrum of the density field is in- 
deed nearly Kolmogorov as well. Another point is that the 
power-spectrum of the density in the isothermal simula- 
tion performed in Sect. 5.2 (let us remember that this sim- 
ulation has exactly the same initial, boundary and forcing 
conditions than the other runs performed with cooling) is 
indeed very close to the Kolmogorov spectrum. The rea- 
son is that the isothermal simulation is subsonic (as for the 
simulation with cooling, the rms velocity is about half the 
sound speed of the WNM). Since in the ISM, the WNM is 
likely to be more or less isothermal until the thermal pres- 
sure is triggered to a value for which WNM is thermally 
unstable, this should constitute a reasonable description 
of WNM before the thermal transition occurs. Therefore, 
we assume that S 2 (k) cx k~ n , with n ~ 8/3. 

To calculate the width of the distribution of the density 
fluctuations for the case of the turbulent WNM, we choose 
the simple window function sharply truncated in k-space, 
W k {R) = OiR- 1 - k) where d(z) = 1 if z > and 
otherwise. With Eq. ® it is easy to see that the integral 
diverges for small k. This is a natural consequence of the 
well known property of the Kolmogorov spectrum that 
most of the energy is on the large scales. We therefore have 
to introduce a lower value in k-spacc, which corresponds 
to an upper limit of the spatial scales, L c . In the case of 
the simulation, this is naturally given by the size of the 
computational box whereas in the ISM, this would be the 
energy injection scale. We have: 

^( fl) =£" Ct ™ = ^((i-)-"-(I)-")(7, 

where C is a constant. 

The question is then, are the density fluctuations rea- 
sonably described by Eqs.[5]and[7]? To answer this ques- 
tion, we have smoothed the density field of the isothermal 
and cooling simulations using the simple window func- 
tion, W b {pi,j) = 1/b 2 T,\i- m \<b,\j-n\<b Pm,n, for various 
values of b. Figures [19] and [20] show the results in the 
isothermal case and in the case with cooling (since the 
isothermal simulation has 5000 2 cells, the values of b are 
2 times lower than in the case with cooling). It is seen 
that, in the isothermal case, the probability distribution 
is reasonably Gaussian. As expected from Eq. [7] a ap- 
pears to be a decreasing function of b = R/dx. In the case 
with cooling, the situation is slightly less clear because the 
shape is a little skewed towards the large densities. This 
is a consequence of the fact that a fraction of the gas is 
thermally unstable and in transition towards CNM or is 
already CNM. This dense gas is therefore on the process 
of condensation or already condensed and should not be 
considered in the pdf of density fluctuations of WNM. 

Beyond this qualitative agreement, we present in 
Appendix A a more quantitative estimate that confirms 
the validity of Eq. [7] Thus it appears that the assumption 
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of the distribution of the density fluctuations being given 
by Eqs.[5]and [7] is sufficiently justified. 

7.2. Mass function 

To calculate the mass function, we can now follow the 
same procedure as for the cosmological calculation. We 
note, however, two important differences. First, since we 
assume statistical stationarity, there is no time depen- 
dence in our analysis. Second, the collapse is not due to 
gravity but to thermal instability. 

Here, since our primary goal is to understand the sim- 
ulation results, we consider the 2D case. Thus, for a per- 
turbation of size, L, the final mass of the structures is 
about M(L) = pL 2 , Eq. [7] is equivalent to 

(7 2 (M) = -^(Af(i c )(«- 2 )/ 2 - M (n-2)/2j-(2-n)/2 ) (g) 

where M = M(R). 

The fraction of objects having a mass larger than M 
is thus given by: 



F(M) 



P R (S)s(5)dS, 



(9) 



where s is the selection function 0. The selection function, 
s, determines the fluctuations that will eventually lead 
to the formation of a CNM structure. For simplicity, we 
will assume that once the density reaches a threshold 5 C , it 
leads to a thermal collapse and to the formation of a CNM 
structure. Although this assumption is reasonable given 
the physical characteristics of thermal instability, we note 
that various improvements could certainly be possible, as 
for example taking into account the stabilizing role of the 
turbulence within the perturbation (in that case s would 
also depend on the velocity dispersion). Thus we have 



/>oo 

F(M) = J P R 



(6)dS. 



The mass function, M(M). 
-dF/dM/{M/p). We have El 



(10) 

is then simply given by 



Af(M)dM = -■ 



1 



P 



S r da 



f^Ma 2 dM 
With Eq. m we get 



exp 



2a 2 



dM. 



da 
dM 



— = -VCp 



-l/2-n/4 



(11) 



(12) 



Af(™- 2 )/ 4 

2M 



(n 



(n-2)/2 _ M (n-2)/2 



2 We note that in cosmology this expression is known to be 
inconsistent because in the cosmological case, the integral of 
dF = f(m)dm from to oo should be equal to 1. This prob- 
lem has been solved by the Excursion set theory (Bond et al. 
1991). The results for the shape of the mass spectrum is how- 
ever independent of this and we will not consider this problem 
further at this stage. 

3 The Excursion set theory would give the same result mul- 
tiplied by a factor 2. 



This leads to: 

N(M)dM <x M {n - 2 ^ 2 - 2 \cx V {-5 2 c l{2a 2 ))dM (13) 

a 

Therefore for small masses, i.e. mass smaller than M(L C ), 
we obtain that Af(M) cx Af(™- 2 )/ 2 - 2 with n ~ 8/3, we 
get Af(M) cx M -5 / 3 . We note that 5/3 ~ 1.666 is rather 
close to the value 1.7 obtained in the simulation. For a 
mass close to M(L C ), there is an exponential cutoff which 
is due to the cut, L c introduced at large scales. Altogether 
these results are similar to the results obtained in the cos- 
mological case. An important difference however lies in the 
exponent of the mass spectrum which is less stiff. This is 
a consequence of the fact that in a turbulent flow, most 
of the energy is at large scales implying more large scale 
fluctuations. 

7.3. The 3D case: similarity with the CO clumps 

This analysis can be straightforwardly applied to the 3D 
case. The only differences are that n = 11/3 (since the 
simulations presented here are bidimensional, it is worth 
to remember that Kim & Ryu (2005) do find that in tran- 
sonic isothermal simulations, the density power spectrum 
is nearly Kolmogorov) and M = ~pR 3 . We obtain: 



da M(™- 3 )/ 6 
dM * 3A~T 



(n-3) 



\ (m(L c )(™" 3 )/ 3 - Af (™- 3 )/ 3 ) 



(14) 



and 



M(M)dM cx Af ( "- 3)/3 ^ 2 — exp(-<5 2 /(2cr 2 ))dAf. (15) 

a A 

The index of the power law of the mass spectrum in the 
3D case and for masses small compared to M(L C ), is about 
(n — 3)/3 ~ 1.78, i.e. slightly stiffer than the value inferred 
in the 2D case. 

Interestingly, we note that the value inferred here, 
turns out to be similar to the index of the mass spec- 
trum observed for the CO clumps of masses between 10 4 
solar mass and one Jupiter mass (Kramer et al. 1998, 
Heithausen et al. 1998). 

Although it cannot be excluded that this is a pure coin- 
cidence, one possibility is that the origin of the CO clumps 
is indeed rooted in the atomic gas. The mass distribution 
of the CO clumps may therefore be determined during the 
atomic phase before the gas becomes molecular. 

8. Conclusion 

We have presented high resolution numerical simulations 
aiming to describe a turbulent interstellar atomic flow. 
The high resolution (2 x 10~ 3 pc) provides a good de- 
scription (although not sufficient to reach numerical con- 
vergence) of the flow close to the spatial scales of the small- 
est structures observed in HI. We confirm the main results 
obtained in paper I. The flow is very fragmented and the 
phases are tightly interwoven. The CNM structures are 
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connected to the WNM by stiff thermal fronts and are 
locally in near pressure equilibrium with the surrounding 
WNM. Altogether, this is reminiscent of the classical 2- 
phase model in spite of the fact that the flow is dynamical 
and turbulent. 

We find that small scale structures, either diffuse or 
very dense, are naturally produced in a turbulent 2-phase 
medium. Whereas the former are simply the tail of the 
structure mass spectrum induced by turbulence, the lat- 
ter, which are somehow reminiscent of the TSAS observed 
in the atomic gas, are transient events produced by super- 
sonic collisions between CNM fragments. 

We characterize the flow by studying various statisti- 
cal quantities, namely the mass distribution of the CNM 
structures, the velocity and density power-spectrum and 
the energy spectrum, paying special attention to the nu- 
merical resolution. We stress that the present problem re- 
quires a lot of numerical resolutions to obtain reliable re- 
sults. Typically at least 5000 2 to 10000 2 cells are needed 
and even so, no strict convergence is achieved. 

Our main results are as follows. 

For structures of size larger than about 10 grid cells, 
for which numerical diffusivity is not too important, we 
find a mass spectrum, M(M) oc M _1 . In order to ex- 
plain this result, we have carried out a calculation based 
on Press & Schecter (1974) formalism. Our theory pre- 
dicts N{M) oc M~ 5 / 3 in 2D and N{M) oc M -16 / 9 in 3D. 
One of the main assumptions is that the CNM clumps are 
due to density fluctuations within WNM, which behaves 
as a sub to transonic nearly isothermal gas. We note that 
these mass spectra are similar to the mass spectrum in- 
ferred observationally for CO clumps. This is compatible 
with the origin of CO clumps being rooted in the very 
diffuse atomic gas although this does not, by any means, 
constitute a proof. 

Like in highly supersonic isothermal flows, the den- 
sity power-spectrum is rather flat. However, the struc- 
tures being pressure-bounded and not bounded by shocks 
as in supersonic isothermal simulations, they are long- 
living objects. Unlike supersonic isothermal simulations, 
the velocity power-spectrum follows the Kolmogorov pre- 
diction and is dominated by its solenoidal part. At scales 
larger than about 1 pc, most of the energy is in the WNM 
whereas at scales smaller than about 1 pc, it is mostly in 
the CNM. Indeed, the energy is nearly equally distributed 
in k-space for scales ranging between about 2 and 2 x 10 -2 
pc where numerical diffusion becomes dominant. This be- 
haviour is due to the density structure of the flow that 
presents relatively flat density power-spectrum, indicating 
that in a thermally bistable flow like the one we studied, 
the kinetic energy at scale below 1 pc, appears to be domi- 
nated by the translational motions of the CNM fragments. 

We also study the influence of the thermal conduction 
by performing various runs with different thermal conduc- 
tivities and find a modest influence, mainly on the small 
structures. 

In a companion paper (paper III) we further character- 
ize the simulation results by studying line of sight statis- 
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Fig. A.l. Isothermal simulation. Density field smoothed 
by the window function Wb — Wib for various value of 
b. The range of the x-axis being proportional to 6™~ 2 , it 
appears that the length of these distribution is reasonably 
gaussian and roughly proportional to b n ~ 2 . 



tics and the physical properties of the CNM clouds. We 
also calculate synthetic HI spectra and reach the conclu- 
sion that, although the lines of sight from which they are 
calculated are very complex, showing in some cases several 
structures, the spectra are relatively smooth and simply 
broadened by the complex motions along the line of sight. 

Altogether these results suggest that the turbulence 
which takes place in the neutral atomic interstellar gas 
is very different from the turbulence which takes place in 
an isothermal or polytropic gas. One important question 
that remains to be investigated is how these conclusions 
will change in a more realistic framework of 3D flows. To 
answer this question 3D simulations have been performed 
and will be presented elsewhere (Audit & Hennebelle 
2007). 
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Appendix A: Further measurement of the pdf of 
density fluctuations 

In order to measure more accurately the dependence of a 
(Eq.[7J on R, we define Sp b - Sp 2 b = W b (p) - W 2 b(p), i.e. 
the density field smoothed by the window function Wb — 
Wzb- With this window function, one finds that a(R) = 
C/(n-2) x (2™~ 2 -l)i?™~ 2 . In particular, it is independent 
of L c and is proportional to R n ~ 2 = i? 2 / 3 for n = 8/3. 

Figures [A. II and lA.21 display the distribution of Spb — 
5p2b for various values of b for the isothermal and the 2- 
phase simulations. The range of the x-axis is proportional 
to b n ~ 2 . As can be seen, the distributions are reasonably 
Gaussian and the width is broadly proportional to b n ~ 2 for 
both cases except for small b in the isothermal simulation. 
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Fig. A. 2. Same as Fig. I A. II for the 2-phase simulation 
in which the density field has been truncated by taking 
min(p, p ). 
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